/**
  CNOK project, Anyang Normal University, IMP-CAS

  \class TAFourier
  \brief This class implements Fourier transfrom of f(r) to f(q), where f(r) is
  spherically symmetric, and q is in x-y plane.
  \author SUN Yazhou, asia.rabbit@163.com
  \since 2020/07/08
  \date Last modified: 2020/07/25 by SUN Yazhou
  \copyright 2020-2023 SUN Yazhou
  \copyright CNOK project, Anyang Normal University, IMP-CAS
*/

#include "config.h"
#ifdef CATCH_CNOK
#include "catch.h"
#endif
#include "TAFourier.h"
#include "TAIntegrate.h"
#include "TAMath.h"
#include "TAException.h"

static const double PI = TAMath::Pi();
static const double FOURPI = 4.*PI;

// computes the Fourier transform of a spherically symmetric function fr,
// defined on a grid of r with nr points. Returns the Fourier  transform in
// function fq on a grid q with nq points.
// the transform is of the form:
// f(q) = \int[dr\vec f(r) exp(-iq\vec.s\vec)] = 4pi/q\int[dr rf(r)sin(qr)]
// integration is over 0->\infty
//
// implements Filon's Integration formula in M Abramowitz and I A Stegun (1970)
// Handbook of Mathematical Functions, Dover, New York, p890, eq.25.4.54
// accuracy~O(h^3) for odd nr (i.e., 0->2n)
// for even nr, the last term is added using trapezoidal rule, with the same
// accuracy (O(h^3))
double TAFourier::Fourier2D(double q, int nr, const double *r, const double *fr){
  // fr is a delta function, so return a normalized const,
  // which is the mass of a single nucleon, i.e. 1, and has been set to nr
  if(!r) return 1.; // or return nr. It's the same.

  if(0. == q){
    double *f = new double[nr]{};
    for(int i = nr; i--;) f[i] = r[i]*r[i]*fr[i]; // r^2*f(r)
    const double fq = FOURPI*TAIntegrate<>::Romberg(nr, r, f);
    delete [] f;
    return fq;
  } // end if
  const double h = (r[nr-1] - r[0])/(nr-1); // r and k step

  const int n = (nr-1)/2; // nr: 0->2n (+1)
  const double k = q; // just a different notation
  const double kh = k*h, kh2 = kh*kh, kh3 = kh2*kh;
  const double sinkh = sin(kh), coskh = cos(kh), sin2kh = 2.*sinkh*coskh;
  const double alpha = 1./kh + sin2kh/(2.*kh2) - 2.*sinkh*sinkh/kh3;
  const double beta  = 2.*((1.+coskh*coskh)/kh2 - sin2kh/kh3); // beta in Filon
  const double gamma = 4.*(sinkh/kh3 - coskh/kh2);
  double sEven = 0, sOdd = 0.; // f(r)sin(kr): 0,2,...,2n; 1,3,..2n-1
  int je, jo; // for even and odd times of h
  for(int j = 1; j <= n; j++){ // j==0 omitted where sin(kr)=0
    je = 2*j; jo = je - 1;
    sEven += je*h*fr[je]*sin(je*kh); // rf(r)sinkr
    sOdd  += jo*h*fr[jo]*sin(jo*kh); // rf(r)sinkr
  } // end for over r
  const double rfrR = je*h*fr[je]; // -- x2n*sin(t2n) --
  sEven -= 0.5*rfrR*sin(je*kh);
  double fq = h*(-alpha*rfrR*cos(kh*2.*n) + beta*sEven + gamma*sOdd);
  // complete the integral with the last interval: Rf(R)sin(kR)*h using trapezoidal rule
  if(0 == nr%2) fq += ( h*(nr-2)*fr[nr-2]*sin(kh*(nr-2))+
    h*(nr-1)*fr[nr-1]*sin(kh*(nr-1)) ) * h / 2.;

  return fq*FOURPI/k;// fq = 4pi/q * \int{dr rf(r)sin(qr)}
} // end member function Fourier2D(...)

/// put the Fourier transform in an array
void TAFourier::Fourier2D(int nr, const double *r, const double *fr, int nq,
    const double *q, double *fq){
  const double h = r[3] - r[2], dk = q[3] - q[2]; // r and k step
  if(h < 0. || dk < 0.)
    TAException::Error("TAFourier", "Fourier2D(...): dr or dq is minus.\
The r and q array must be of ascending order.");

  // calculate the integral for each momentum k
  for(int i = 0; i < nq; i++) fq[i] = Fourier2D(q[i], nr, r, fr);
} // end member function Fourier2D(...)


#ifdef CATCH_CNOK
// UNIT TEST using catch2 //
TEST_CASE("testing Filon's integration formula", "[filon]"){
  const int n = 200;
  double r[n]{}, rho[n]{}, f[n]{}; // f(r) = 4pi/q \int_0^{\pi}{r*rho*sin(qr)}
  const double h = PI / n; // the step
  const double q = 2.;
  // f(x) = sin(qx); I=\int_0^{\pi}{sin(qx)}
  for(int i = 0; i < n; i++){
    r[i] = i*h; rho[i] = 1.;
    f[i] = r[i]*rho[i]*sin(q*r[i]);
  } // end for
  double i1 = TAIntegrate<>::Romberg(n, r, f)*4.*PI/q;
  double i2 = TAFourier::Fourier2D(q, n, r, rho);
  double i3 = TAIntegrate<>::Simpson(n, r, f)*4.*PI/q;
  CHECK(i1 == Approx(-PI*PI).epsilon(1.e-3));
  CHECK(i2 == Approx(-PI*PI).epsilon(1.e-3));
  CHECK(i3 == Approx(-PI*PI).epsilon(1.e-3));
} // end test case
#endif
